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Abstract 

We numerically study the scaling behavior of period doublings at the zero- 
coupling critical point in a four-dimensional volume-preserving map consisting 
of two coupled area-preserving maps. In order to see the fine structure of 
period doublings, we extend the simple one-term scaling law to a two-term 
scaling law. Thus we find a new scaling factor 5^ (= 1.8505...) associated 
with scaling of the coupling parameter, in addition to the previously known 
scaling factors 6i (= —8.7210...) and 62 (= —4.4038...). These numerical 
results confirm the renormalization results reported by Mao and Greene [Phys. 

Rev. A 35, 3911 (1987)]. 

PACS numbers: 05.45.+b, 03.20.-hi, 05.70.Jk 
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Universal scaling behavior of period doubling has been found in area-preserving maps 
||T|-^. As a nonlinearity parameter is varied, an initially stable periodic orbit may lose its 
stability and give rise to the birth of a stable period-doubled orbit. An infinite sequence of 
such bifurcations accumulates at a finite parameter value and exhibits a universal limiting 
behavior. However these limiting scaling behaviors are different from those for the one- 
dimensional dissipative case P]. 

An interesting question is whether the scaling results of area-preserving maps carry over 
higher-dimensional volume-preserving maps. Thus period doubling in four-dimensional (4D) 
volume-preserving maps has been much studied in recent years [|7|,|9Hl3|. It has been found 



in Refs. ||TTHl3[ that the critical scaling behaviors of period doublings for two symmetrically 



coupled area-preserving maps are much richer than those for the uncoupled area-preserving 
case. There exist an infinite number of critical points in the space of the nonlinearity and 



coupling parameters. It has been numerically found in [|ri|)0 that the critical behaviors 
at those critical points are characterized by two scaling factors, 6i and 62- The value of 61 
associated with scaling of the nonlinearity parameter is always the same as that of the scaling 
factor 6 (= 8.721 . . .) for the area-preserving maps. However the values of 62 associated with 
scaling of the coupling parameter vary depending on the type of bifurcation routes to the 
critical points. 



The numerical results [11,12| agree well with an approximate analytic renormalization 



results obtained by Mao and Greene |jT3[ , except for the zero-coupling case in which the two 



area-preserving maps become uncoupled. Using an approximate renormalization method 
including truncation, they found three relevant eigenvalues, 61 = 8.9474, 62 = —4.4510 and 
^3 = 1.8762 for the zero-coupling case |]14[. However they believed that the third one 63 



is an artifact of the truncation, because only two relevant eigenvalues 61 and 62 could be 
indentified with the scaling factors numerically found. 

In this paper we numerically study the critical behavior at the zero-coupling point in 
two symmetrically coupled area-preserving maps and resolve the inconsistency between the 
numerical results on the scaling of the coupling parameter and the approximate renormaliza- 



tion results for the zero-coupling case. In order to see the fine structure of period doublings, 
we extend the simple one-term scaling law to a two-term scaling law. Thus we find a new 
scaling factor ^3 = 1.8505 . . . associated with coupling, in addition to the previously known 
coupling scaling factor 82 = —4.4038 .... The numerical values of 62 and ^3 are close to 
the renormalization results of the relevant coupling eigenvalues 62 and ^3. Consequently 
the fixed map governing the critical behavior at the zero-coupling point has two relevant 
coupling eigenvalues 62 and ^3 associated with coupling perturbations, unlike the cases of 
other critical points. 

Consider a 4D volume-preserving map consisting of two symmetrically coupled area- 
preserving Henon maps with a periodic boundary condition. 



T : <^ 



x(t + l) = -y(t) + F(x(t)), 
y(t+l) = x(t). 



where x = (xi,X2), y = (?/i,l/2), F(x) = (Fi(x), F2(x)). Here z„(t) = (^^(t), ^^(t)) 
(m = 1,2) is the state vector of the mth element at a discrete time t, and the periodic 
boundary condition imposes Zm = Zm+2 for all m. The mth component -Fm(x) of the vector- 
valued function F(x) is given by 

Fm{x) = F{Xm, Xm+l) = f{Xm) + g{Xm, Xm+l), (2) 

where / is the nonlinear function of the uncoupled Henon's quadratic map [|l^, i.e., 

f{x) = l-ax\ (3) 
and g{xi, X2) is a coupling function obeying a condition 

g{x,x) = for any x. (4) 



The two-coupled map (|l]) is called a symmetric map [|rT|,0 because it has an exchange 
symmetry such that 

a-^Ta{z) = T(z) for all z, (5) 



where z = {zi, Z2), o is an exchange operator acting on z such that crz = (2:2, zi), and 
is its inverse. The set of all fixed points of a is a two-dimensional (2D) subspace of the 
4D state space, on which x\ = X2 and yi = y2- It follows from Eq. (^) that the exchange 
operator a commutes with the symmetric map T, i.e., aT = Ta. Hence the 2D subspace 
becomes invariant under T. An orbit is called an in-phase orbit if it lies on the 2D invariant 
subspace, i.e., it satisfies 



Xi{t) = X2{t) = x{t), yi{t) = y2{t) = y{t) for all t. 



(6) 



Otherwise it is called an out-of-phase orbit. Here we study only in-phase orbits. They can 
be easily found from the uncoupled Henon map since the coupling function g obeys the 
condition @. 

The Jacobian matrix DT of the two-coupled map T is: 

r\ 



DT 



DF 
I 







(7) 



where DF is the Jacobian matrix of the function F(x), / is the 2x2 identity matrix, and 
is the 2x2 null matrix. Since Det{DT) = 1, the map T is a 4D volume-preserving map. 
Furthermore, if DF is a symmetric matrix, i.e., DF* = DF (t denotes transpose), then the 
map T is a symplectic map because its Jacobian matrix satisfies the relation DT^ J DT = J 
TB| , where 

\ 

(8) 



J 



/ 

V 



/ 
-/ 

Stability analysis of an orbit in coupled maps is conveniently carried out by Fourier- 
transforming with respect to the discrete space {m} |T^. Consider an orbit {x(t)} = 
{xm(t) ; m = 1, 2} of the two-coupled map T expressed in the form of second-order difference 
equations. 



T : x^{t + 1) + x^t - 1) = FUm) 

= f{Xm{t)) + g{xm{t),Xm+i(t)), m = 1, 2. 



(9) 



The discrete spatial Fourier transform of the orbit is: 



T[xm{t)] = -^e--^™^x^(t) = J = 0, 1. 

^m=l 



(10) 



The wavelength of a mode with index J is 

To determine the stability of an in-phase orbit [xi{t) = X2{t) = x(t) for all t], we consider 
an infinitesimal perturbation {Sxm{t)} to the in-phase orbit, i.e., Xm{t) = x{t) + 6xm{t) for 
m = 1,2. Linearizing the two-coupled map (P) at the in-phase orbit, we have 



N 



5Xm{t + 1) + 5Xm{t - 1) = f{x{t)) 5Xm{t) + Y^Gi{x{t)) 



(11) 



where 

' X\=X2 = X 

Hereafter we will call the functions Gis "reduced" coupling functions of ^'(x) 
Let 5^j{t) be the Fourier transform of dxmit), i-e.. 



(12) 



S^t) = n^x^it)] = 7;i: e-"'"Mx^(t), J = 0, 1. 

^m=l 



Then the Fourier transform of Eq. ([TlD becomes: 



(13) 



5Ut + 1) + SUt - 1) = [f'{x{t)) + Y.Giix{t))e 

1=1 



m{l-l)j] 



5^,{t), J = 0,1. 



This equation can be also put into the following form: 



Hi) 



J 



\6^,{t-l) ) 



J = 0,1, 



where 



L,(t) 



h'{x{t)) + Y^Gi{x{t))e-^('-^~^^ 



1=1 



\ 







Note that the determinant of Lj is one, i.e., Det{Lj) = 1. 

It follows from the condition (H) that the reduced coupling functions satisfy 



(14) 



(15) 



(16) 



T.Gi{x) = 0. (17) 

1=1 

Hence there exists only one independent reduced coupling function G{x) such that 

G2{x) = G{x), Gi{x) = -G{x). (18) 



Substituting Gi's into the (1, 1) entry of the matrix Lj{t), we have: 

^ n for j = 0, 

Y^Giixit))e''^^'''^^ 

-2G{x{t)) for J = 1 



(19) 



Hence the matrices L^-'s of Eq. ([Iq ) are real matrices. 

Stability of an in-phase orbit of period q is determined by iterating Eq. ([T5|) q times: 
/ 



6Ut + q) 

S^j{t + q-l) 



where 



= M, 












/ 








t+q- 


-1 







J = 0,1, 



(20) 



M, = n L,ik)- 

k=t 



(21) 



That is, the stability of each mode with index j is determined by the 2x2 matrix Mj. 
Since Det{Mj) = 1, each matrix Mj has a reciprocal pair of eigenvalues, \j and \J^. These 
eigenvalues are called the stability multipliers of the mode with index j. Associate with a 
pair of multipliers (Aj, A~^) a stability index, 



Pj = -^i + \ J = 0, 1 



(22) 



which is just the trace of Mj, i.e., pj = Tr{Mj). Since Mj is a real matrix, pj is always real. 
Note that the matrix Mq for the case of the j = mode is just the Jacobian matrix of the 
uncoupled Henon map. Hence the coupling affects only the stability index pi of the j = 1 
mode. 

It follows from the reality of pj that the reciprocal pair of eigenvalues of Mj lies either on 
the unit circle, or on the real line in the complex plane, i.e., they are a complex conjugate 
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pair on the unit circle, or a reciprocal pair of reals. Each mode with index j is stable if 
and only if its stability index pj is real with \pj\ < 2, i.e., its stability multipliers are a pair 
of complex conjugate numbers of modulus unity. A period-doubling (tangent) bifurcation 
occurs when the stability index pj decreases (increases) through —2 (2), i.e., two eigenvalues 
coalesce at Xj = —1 (1) and split along the negative (positive) real axis. 

An in-phase orbit is stable only when all its modes are stable. Hence its stable region in 
the space of the nonlinearity and coupling parameters is bounded by four bifurcation lines 
associated with tangent and period-doubling bifurcations of both modes (i.e., those curves 
determined by the equations pj = ±2 for j = 0, 1). When the stability index po decreases 
through —2, the in-phase orbit loses its stability via in-phase period-doubling bifurcation 
and gives rise to the birth of the period-doubled in-phase orbit. Here we are interested in 
scaling behaviors of such in-phase period-doubling bifurcations. 

As an example we consider a linearly-coupled case in which the coupling function is 

g{xi,X2) = -{x2- xi). (23) 

Here c is a coupling parameter. Figure |T] shows the stability diagram of in-phase orbits 
with period g = 1,2, and 4 [|T^]. As previously observed in [^,0, each "mother" stability 



region bifurcates into two "daughter" stability regions successively in the parameter plane; 
hereafter we call the direction of the left (right) branch of the two daughter stability regions 
L (R) direction. Consequently the stability region of the in-phase orbit with period q = 2^ 
{n = 0,1,2,.. .) consists of 2" branches. Each branch can be represented by its address 
[oq, . . . , a„], which is a sequence of symbols L and R such that = L and ai = L or R for 
i > 1. 

An infinite sequence of connected stablity branches (with increasing period) is called a 



bifurcation "route" |Tl],|T2[ . Each bifurcation route is also represented by an infinite sequence 
of symbols L and R. A "self-similar" bifurcation "path" in a bifurcation route is formed by 
following a sequence of parameters (a„, c„), at which the in-phase orbit of level n (period 2"^) 
has some given stability indices (po)Pi) (e.g., Po = —2 and pi = 2) [|Tli|l2| . All bifurcation 



paths within a bifurcation route converge to an accumulation point (a*, c*), where the value 
of a* is always the same as that of the accumulation point for the area-preserving case (i.e., 
a* = 4.136166 803 904. . .), but the value of c* varies depending on the bifurcation routes. 
Thus each bifurcation route ends at a critical point (a*,c*) in the parameter plane. 

It has been numerically found that scaling behaviors near a critical point are characterized 
by two scaling factors, 6i and 62 [TllJl^. The value of 61 associated with scaling of the 



nonlinearity parameter is always the same as that of the scaling factor 6 (= 8.721 .. .) for 
the area-preserving case. However the values of 62 associated with scaling of the coupling 
parameter vary depending on the type of bifurcation routes. These numerical results agree 



well with analytic renormalization results [13|, except for the case of one specific bifurcation 



route, called the E route. The address of the E route is [(L, R, )°°] (= [L, R, L,R, . . .]) and 
it ends at the zero-coupling critical point (a*,0). 

Using an approximate renormalization method including truncation, Mao and Greene 



Ij] obtained three relevant eigenvalues, 61 = 8.9474, 62 = —4.4510, and ^3 = 1.8762 for 
the zero-coupling case; hereafter the two eigenvalues 62 and ^3 associated with coupling 
will be called the coupling eigenvalues (CE's). The two eigenvalues 61 and 62 are close to 
the numerical results of the nonlinearity-parameter scaling factor 6i{= 8.721 . . .) and the 
coupling-parameter scaling factor 62{= —4.403 . . .) for the E route. However they believed 
that the second relevant CE ^3 is an artifact of the truncation, because it could not be 
identified with anything obtained by a direct numerical method. 

In order to resolve the inconsistency between the numerical results and the renormal- 
ization results for the zero-coupling case, we numerically reexamine the scaling behavior 
associated with coupling. Extending the simple one-term scaling law to a two-term scaling 
law, we find a new scaling factor ^3 = 1.8505 . . . associated with coupling in addition to the 
previously found coupling scaling factor 62 = —4.4038 . . ., as will be seen below. The values 
of these two coupling scaling factors are close to the renormalization results of the relevant 
CE's 62 and ^3. 

We follow the in-phase orbits of period 2" up to level n = 14 in the E route and 
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obtain a self-similar sequence of parameters (an,c„), at which the pair of stability indices, 

{Po,n,Pi,n), of the orbit of level n is (—2,2). The scalar sequences {a„} and {c„} converge 

geometrically to their limit values, a* and 0, respectively. In order to see their convergence, 

define 5„ = ""'"'^ and /i„ = ""^"^ , where Aa^ = — a„-i and Ac„ = Cn — c„_i. Then 
Aa„ Ac„ 

they converge to their limit values 6 and /i as n — > oo, respectively. Hence the two sequences 
{Aa„} and {Ac„} obey one-term scaling laws asymptotically: 

Aa„ = C(")r", Ac„ = C('=V~" for large n, (24) 

where C*-"^ and C*^^-* are some constants. 

The two sequences {6n} and {pn} are shown in Table |. The second column shows rapid 
convergence of 5„ to its limit value S (= 8.721 . . .), which is close to the renormalization 
result of the first relevant eigenvalue (i.e., 6i = 8.9474). The sequence {/in} also seems to 
converge to fi = —4.403 . . . (see the third column). However its convergence is not as fast 
as that for the case of the sequence {5„}. This is because two relevant CE's 62 and ^3 are 
involved in the scaling of the sequence 

Taking into account the effect of the second relevant CE ^3 on the scaling of the sequence 
{Ac„}, we extend the simple one-term scaling law (plj) to a two-term scaling law: 

Acn = Cifi]^" + 6*2/12" for large n, (25) 

where |/ii| > |/i2|- This is a kind of multiple scaling law |T^,2^. Eq. (^) gives 



Ac„ — tiAc„_|_i — /:2Ac„_|_2, (26) 

where ti = /ii + /i2 and ^2 = ^1^2- Then /ii and /i2 are solutions of the following quadratic 
equation, 

/i^ - ti/i + t2 = 0. (27) 
To evaluate fii and /i2, we first obtain ti and ^2 from Aa„'s using Eq. (|26D: 

^ _ ACnACn+l - ACn-iACn+2 ^ _ Ac^ - ACn+lACn-l .^^^ 

^ A4^i - Ac„ Ac„+2 ' ^ Acl^^ - Ac„ Ac„+2 ' 



Note that Eqs. (p5D-(p8D hold only for large n. In fact the values of tj's and /ij's {i = 1,2) 
depend on the level n. Therefore we explicitly denote tj's and /Xj's by tj,n's and iii^n^i 
respectively. Then each of them converges to a constant as n — > oo: 

lim ti n = ti, lim /Xj^ = /Xj, i = 1, 2. (29) 

n^oo ' n—>-oo ' ' ' 

Three sequences {/Ui,„}, {fJ'2,n}, and {fil J /^2,n} are shown in Table 0. The second column 
shows rapid convergence of /ii_„, to its limit values /ii (= —4.403 897805). Comparing this 
sequence with the sequence {/in} in Table I, we find that the finite sequence {Ac„} obeys 
the two-term scaling law better than the one-term scaling law. That is, the accuracy for 
the scaling factor /ii obtained from the two-term scaling law is much higher than that 
obtained from the one-term scaling law. The value of this scaling factor /ii is close to the 
renormalization result of the first relevant CE (i.e., = —4.4510). From the third and 
fourth columns, we also find that the second scaling factor fi2 is given by a product of two 
relevant CE's 62 and ^3, 

/i2 = |, (30) 

where 62 = /ii and ^3 = 1.850 65 . It has been known that every scaling factor in the multiple- 
scaling expansion of a parameter is expressed by a product of the eigenvalues of a linearized 
renormalization operator |T9|j20|| . Note that the value of ^3 is close to the renormalization 



result of the second relevant CE (i.e., ^3 = 1.8762). 

We now study the coupling effect on the stability index pi „ of the j = I mode of the 
in-phase orbit of period 2" near the zero-coupling critical point (a*, 0). Figure ^ shows three 
plots of pi,„(a*,c) versus c for n = 4,5, and 6. For c = 0, pi,„ converges to a constant pi 



-2.543 510 20. . .), called the critical stability index ||12[, as n ^ 00. However, when c 

^^'"^ ) at the zero-coupling 

{a*,0) 



is non-zero pi „ diverges as n — > cxd, i.e., its slope Sn (= „ 

oc 

critical point diverges as n 00. 



The sequence {Sn} obeys a two-term scaling law, 

Sn = Diiy^ + D2iy^ for large n, (31) 
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where > |z/2|. This equation gives 

Sn+2 = T^lSn+l — ^28^ (32) 

where ri = ui + 1/2 and r2 = i'ii'2- As in the scahng for the couphng parameter, we first 
obtain ri and r2 of level n from S'„'s: 

C C C C C2 C C 

~ Q2 C C ' ^2,n — ^; • \^^) 

'-'n.+l'-'n-l 

Then the scahng factors z/i^„ and z/2,„ of level n are given by the roots of the quadratic 
equation, — ri „z/„ + r2^„ = 0. They are listed in Table |T| and converge to constants ui 
(= —4.403 897805 09) and z/2 (= 1.850 535) as n ^ 00, whose accuracies are higher than 
those of the coupling-parameter scaling factors. Note that the values of ui and z/2 are also 
close to the renormalization results of the two relevant CE's 62 and 53. 

We have also studied several other coupling cases with the coupling function, g{xi,X2) = 
-(^2 — x") (n is a positive integer). In all cases studied {n = 2,3,4,5), the scaling factors 
of both the coupling parameter c and the slope of the stability index pi are found to be the 
same as those for the above linearly-coupled case {n = 1) within numerical accuracy. Hence 
universality also seems to be well obeyed. 
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TABLES 

TABLE I. Scaling factors 6n and /x„ in the one-term scaling for the nonlinearity and coupling 

parameters are listed. 



n 






6 


8.721086 300 39 


-4.399 51 


7 


8.721096 265 95 


-4.405 75 


8 


8.721097056 68 


-4.403 12 


9 


8.721 097 188 39 


-4.404 22 


10 


8.721097198 70 


-4.403 76 


11 


8.721097200 44 


-4.403 96 


12 


8.721097200 58 


-4.403 87 


13 


8.721097 200 60 


-4.403 91 



TABLE II. Scaling factors and fi2,n in the two-term scaling for the coupling parameter 



are shown in the second and third columns, respectively. A product of them, — '—, is shown in the 
fourth column. 



n 




M2,n 


Ml,n 
^2,n 


5 


-4.403 908 128 


10.4374 


1.85817 


6 


-4.403 899 694 


10.465 9 


1.853 09 


7 


-4.403 898 736 


10.458 2 


1.85446 


8 


-4.403 897 867 


10.474 8 


1.85152 


9 


-4.403 897847 


10.473 9 


1.85168 


10 


-4.403 897806 


10.478 4 


1.850 89 


11 


-4.403 897 807 


10.478 6 


1.850 85 


12 


-4.403 897805 


10.479 7 


1.850 65 
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TABLE III. 


Scaling 


factors vi^n and i'2,n in the two-term scaling 


; for the slope of the stability 


indc}x (^f tli(; j = 


1 mode; 






n 








5 




-4.403 898453 59 


1.851433 5 


6 




-4.403 897730 29 


1.850 782 6 


7 




-4.403 897813 85 


1.850 603 6 


8 




-4.403 897804 07 


1.850 553 8 


9 




-4.403 897805 21 


1.850 540 


10 




-4.403 897805 07 


1.850 5361 


11 




-4.403 897805 09 


1.850 535 


12 




-4.403 897805 09 


1.850 534 9 
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FIGURES 

FIG. 1. Stability diagram of in-phase orbits for the linearly-coupled case. The horizontal 

(non-horizontal) solid and dashed lines denote the period-doubling and tangent bifurcation lines 
of the j = (1) mode, respectively. For other details see the text. 



FIG. 2. Plots of the stability index pi^„(a*, c) versus c for n = 4, 5, 6. 
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